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Abstract 

A reliable and fast numerical scheme is crucial for the ID simulation 
of blood flow. In this paper, a ID blood flow model is incorporated with 
a Kelvin- Voigt viscoelastic constitutive relation of wall. This lead to a 
nonlinear hyperbolic-parabolic system, which is then solved with four nu- 
merical schemes, namely: MacCormack, Taylor-Galerkin, second order 
finite volume and local discontinuous Galerkin. The numerical schemes 
are tested on an uniform vessel, a simple bifurcation and a network with 
55 arteries. In all of the cases, the numerical solutions are checked fa- 
vorably against analytic, semi-analytic solutions or clinical observations. 
Among the numerical schemes, comparisons are made in four aspects: the 
accuracy, the ability to capture shock-like phenomena, the computation 
speed and the complexity of the implementation. The suitable conditions 
for the application of each scheme are discussed. 

Keywords: blood flow; ID flow modeling; vascular network; numerical simu- 
lation 



1 Introduction 

There exist direct 3D simulations of blood flow in arteries, nevertheless they 
are known to be time and memory consuming and therefore most of them are 
restricted to local positions (single vessel, confluences) [HE]. If we assume an 
axisymmetric circular velocity profile in the vessel, the 3D problem can be re- 
duced to a 2D problem. If we further take the advantage of long wave length, a 
ID model can be obtained. The ID model is specially interesting for several rea- 
sons. First, this model captures well the behaviours of pulse wave, from which 
one can extract a lot of useful information about the cardiovascular system. For 
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example, the Pulse Wave Velocity (PWV) has been recognized by European 
Society of Hypertension as a very important marker to the diagnosis and treat- 
ment of hypertension |24[ [5j. Second, it allows fast numerical computations, 
which allows real-time applications for medical planning. Third, it also provide 
pertinent boundary conditions for 3D simulations in multi-scale models [S]. 

The ID governing equations can be obtained by integrating Navier-Stokes 
equations across the cross-section of a circular vessel with the assumption of a 
long- wave perturbation and axisymmetric velocity profile [HI [38l Q~B] ■ It results 
in a system of two partial differential equations (PDEs) for mass and momentum 
conservation which involve the flow rate Q, the cross-section area A and the 
average pressure P. To close the system, the constitutive relation of the arterial 
wall comes in, which relates P and A. After the insertion of this relation into 
the PDEs, a nonlinear hyperbolicity-dominated system is obtained. Depending 
on the details of the model of the arterial wall, there may be a diffusive term 
by the viscosity of the wall or/and a dispersive term by the axial tension. 

In case of weak nonlinearity (i.e. small perturbation around the equilibrium 
state [H1[3I]), we can linearize the ID governing equations and find solutions in 
frequency domain [29, 42 . But for the full nonlinear system, analytic solutions 
are not available yet. Thus several numerical schemes have been proposed and 
used to solve the system in the time domain, and roughly we classify them in: 

• Finite Difference (FD) [36l E3 E21 GH HQ1 SB] 

• Finite Volume (FV) [7J El HE] 

• Finite Element (FE), [23 [Q [371 US] 

• Discontinuous Galerkin (DG) [H [26j [28l EH [37] 

These schemes have been successfully applied in other communities where 
researchers have to solve similar hyperbolic problems. For instance, the Mac- 
Comack scheme (FD) was principally designed for gas dynamics (i.e. ID com- 
pressible Euler equations) and it was then successfully used to compute blood 
flow in veins |12j . From ideas frequently applied in shallow water equations, 
Delestre et al. [7J obtained "well balanced" schemes which properly treat the 
source term induced by a tapered artery. The ID model and the numerical solu- 
tions have been validated by in vitro experimental [36, 43, 1 or in vivo clinical 
data [33] [M] . But usually only one particular scheme was chosen in a study 
and no cross comparisons among the schemes can be found. Thus the advan- 
tage/drawback and the suitable conditions for applications of each scheme were 
not discussed in those literatures. 

Our objective in this paper is to make a cross comparison of the four kind of 
schemes and to suggest the suitable conditions of application for each scheme. 
In general, we note that FD schemes are not flexible to treat complex compu- 
tational geometries in high dimensions (2D or 3D). However, in ID dimension 
and low order accuracy schemes, FD, FE and FV schemes are in fact completely 
equivalent for linear problems. But for problems with large nonlinearities, solu- 
tions with sharp gradient may appear and the performance of different schemes 
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could be different. Equally important is the accuracy. For DG scheme it may 
be tuned either by the degree of the polynomial or by the mesh size. But if 
a diffusive term is added to the governing equations, it will be hard to treat 
by an implicit time marching method (e.g. Crank-Nicolson) in the DG setting, 
thus the time step may be very severely limited for stability. Therefore, the 
performance of each scheme depend on the main features of the problems at 
hand. And the problems with different main features actually arise in a wide 
range of applications. For instance, no shock is observed in arteries in normal 
physiological conditions but shock-like phenomena may arise in veins [25] or 
in arteries when the human body suffers from a blunt impact by accident |16j . 
For another instance, in some conditions diffusive term or dispersive term may 
arise as a source term [T] and the proper treatment of these terms will pose 
different levels of difficulty in each numerical framework. Thus to make a cross 
comparison of the numerical schemes is interesting and useful. 

In this paper, Section [2] presents the governing equations and the charac- 
teristic structure of the homogeneous part of the nonlinear system. Section [3] 
presents the numerical solvers. In particular, a large amount of details of com- 
putation are given because this kind of information is scattered in literature. 
In this section, firstly an operator splitting is proposed (in the FD, FV and FE 
frameworks) to separate the hyperbolic and parabolic part. Then the treat- 
ment of the boundary conditions is discussed. Following that, MacCormack, 
Taylor-Galerkin, and second order finite volume are presented to integrate the 
hyperbolic subproblem. The parabolic subproblem is treated by Crank-Nicolson 
method. At the end of this section, a local discontinuous Galerkin method is 
presented for the hyperbolic-parabolic problem. Section [2] shows the analytic 
and numerical results. The system is linearized and asymptotic solutions are 
obtained for cases when there are different source terms in the system. The 
physical interpretations of the solutions are various behaviours of the wave: 
propagation, attenuation or diffusion in an uniform tube. The effect of the skin 
friction and the viscosity of the wall on the pulse wave is clearly observed. More- 
over, a case with a larger nonlinearity is computed and the ability of the four 
schemes to properly capture the shock-like phenomena is tested. After the tests 
on a single vessel, a simple bifurcation is computed and the numerical reflection 
and transmission coefficients are compared with analytic ones predicted with 
linearized equations. Finally, a network with 55 arteries is computed. In all 
of the cases, the numerical solutions are compared favorably with the analytic, 
semi-analytic solutions or clinical observations. In the last section, comparisons 
among the four schemes are made in four aspects: the accuracy, the ability to 
capture shock-like phenomena, the computation speed and the complexity of 
the implementation. The suitable conditions for the application of each scheme 
are discussed. 
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2 The ID model of arterial blood flow 



2.1 ID mathematical model 

Under the assumption of an axisymmetric velocity profile, the ID arterial blood 
flow model can be written as: 

^ + ^ = °< (la) 
dQ ^ d Q 2 .AdP Q 
^ + d- x {a ^ ) + ^ = - Cf A> (lb) 

where as stated above, A is the cross-section area of the artery, Q the flow 
rate or flux, and p the density of blood. The coefficient a is the momentum 
correction factor, and Cf is the skin friction coefficient. They depend on the 
shape of the velocity profile. Usually, the profile can be estimated from the 
Womersley number which is defined as R^uj/u, with R the radius of the vessel, 
ui the frequency of the pulse wave and v the kinematic viscosity of the fluid. 
With a small Womersley number, we can take a Poiseuille (parabolic) profile. 
In that case a = | and Cf = 8ttv. This choice is only valid for very viscous 
flows [HIUS]. In practice, viscosity is not so large, and the profile is more flat. 
For a completely flat profile a equals 1. This value is often used since it leads to 
a considerable simplification in analysis and the loss of relevance of the model is 
very small in most cases 10J . Thus we assume its value is I in this paper. The 
value of Cf needs special attention because it has significant influence on the 
pulse wave. In practical applications, its value has to be determined according 
to the particular problem at hand (both in vitro and in vivo). We assume its 
value is 8ttv according to a Poiseuille profile. We are aware of the limit of this 
approximation. However, as our purpose is comparison of numerical schemes, 
we do not discuss any more the value of a and Cf. 

To close the system, several viscoelastic constitutive relations for arterial 
wall have been presented in literature, like [TJ [33]. We choose the Kelvin- 
Voigt model for simplicity. We assume that the arterial wall is thin, isotropic, 
homogeneous, incompressible, and moreover that it deforms axisymmetrically 
with each circular cross-section independently of the others. We denote the 
undeformed cross-section area by A and the external pressure of the vessel by 
P ex t- Then, the relation linking A and P is: 

dA 



P = P ext +f3(VA-^A ) + v s — , (2) 
with the stiffness coefficient (3, 

\FkE\i 



and the viscosity coefficient v s 



P (i-^W 



2(1-7,2)^4' 



(3) 
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where 77 is the Poisson ratio, which is 0.5 for an incompressible material, E the 
Young's modulus, h the thickness of the wall and 4> the viscosity of the material. 
For convenience, we further define C v = ^f- for reasons which will be clear very 
soon in the next section. We also note that in absence of the wall viscosity we 
retrieve the classical Hooke's law. 



2.2 Characteristic structure of the system 

After presenting the system of equations, we remind its hyperbolic feature by 
discussing the characteristic structure. The discussion is classical, and can be 
found in text books [HJ |2Qj . The notations we introduce here will be useful for 
the discussion of the numerical solvers. We assume P ex t is constant along the 
axial variable x, and substitute the constitutive relation ([2| into Eq. (lb). We 



note that ^ can be replaced by — ^ thanks to Eq. la The equation for the 
balance of momentum turns out 

HQ. + — (91 + P-Ai) - — — ( ggw r — A ( 9 (P^) l^A^h 
dt dx A ip ' pdx {Vs dx } f A dx 3 dx'' 

(4) 

Under the assumption of a small perturbation of A, we approximate the term 
jic(^) h y C v& the already defined coefficient C v = *f = 2p(1 ^^ 
which turns out to be independent of A or Q. The governing equations may be 
written as: 

dU dF S, (5) 



dt dx 



where 



and 







In this equation, U is the conservative variable, F the corresponding flux and S 
the source term. Note that the flux (scaled by constant density) consists of two 
parts, the convective F c and the diffusive F v . We recognize ^- due to the fluid 
flow, J^-A^ due to the elasticity, and —C v ^ due to the viscosity of the wall. In 
generaf, the suitable numerical techniques for the convective and diffusive flux 
are different. Thus it is common to separate the diffusive term and put it on 
the right side. Thus we write now the problem in a convection-diffusion form: 

9U dF „ „ 

lU + 0x- = S + D (6) 



with 



F = F C , D 



dx 2 



5 




We consider firstly the homogeneous part and later the non-homogeneous part. 
Expanding the derivative of the flux, the homogeneous part can be written in a 
quasi-linear form 

dU dU n 
where J c is the Jacobian matrix 

T -( ° 1 

\ A 2 +C Z A 

with the Moens-Korteweg celerity 

(8) 



Actually, A is always positive. Therefore c is real, which is the speed of the 
pressure wave with respect to the fluid flow. The matrix J c has two different 
eigenvalues 

Ai l2 = !±c (9) 

Linear algebra shows J c must be diagonalizable in the form J c = RABT 1 . The 
columns of R are the right eigenvectors of J c . Left multiplying Eq. Q by R~ x , 
one obtains 

at ox 

By introducing a new unknown variable vector which satisfies djjW = R , the 
previous equation can be transformed into 

dW K dW „ 

-m +A l^ = - (10) 

Wtp can be readily obtained by integrating djjW = i? -1 componentwise 

Wi, 2 = Q±4c. (11) 

W = [Wi,W2] T is called Riemann invariant vector or characteristics. In time- 
space plane, are constants along the lines D t Xi^{t) — Xi,2- In physiolog- 
ical conditions, Ai > > A2. The two families of characteristic propagate in 
opposite directions. The homogeneous part is a subcritical hyperbolic system. 



For further use, we get the expression for A and Q by inverting the relation ( 11 ) 



A= {w 1 -W2f M 2 Q = A w 1± w 1 

1024 \P J 2 K ' 

In the non- homogeneous part, the skin friction term dissipates the momen- 
tum and the second order derivative of Q is diffusive. Thus the full system 
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has hyperbolic-parabolic features. In physiological conditions, the Womersley 
number is not too big and the artery is almost uniform, thus the source term 
will be very small and the system is dominated by the hyperbolicity feature. 
If the properties of the artery have sharp variations, large source terms will be 
introduced. In this case, we will treat the artery as two different connected 
segments. 



3 Numerical solvers 

Having defined the problem and notations, in this section we present the nu- 
merical solvers. The original problem is split into two subproblems, hyperbolic 
and parabolic. Three numerical schemes are presented to treat the hyperbolic 
subproblcm. For the parabolic subproblem, Crank-Nicolson method is suitable. 
Because of the duplication of values at the interface of elements in the DG 
setting, there are difficulties to apply Crank-Nicolson scheme. A local discon- 
tinuous Galerkin method is adopted to treat the unsplit problem. 



3.1 Operator splitting 

There are explicit high resolution schemes for hyperbolic problems. But for 
parabolic problem, implicit schemes are necessary in general for a reasonable 
time step for time integration. Thus we applied a fractional step or operator 
splitting method. Starting from Eq. ([6]), the original problem is split into to a 
hyperbolic subproblem, 

and a parabolic one, 

Let us consider the time intervals (t n ,t n+1 ), for n = 0,1,..., with t n = nAt. 
In every time interval, the hyperbolic problem is solved to get a predictor U*, 
which is used as the initial condition (I.C.) of the second problem. The second 
step can be viewed as a corrector. The original problem is approximated by a 
sequential application of the two subproblems in a certain order. 

From data U n , we may make a prediction U* by evolving time At of the 
hyperbolic subproblem, and correct it with the evolution over At of the parabolic 
subproblem, 

AtH AfP 

u n - — ► u* - — ► u n+ \ 

where e Am (e Atv ) means to solve the hyperbolic (parabolic) subproblem over 
At. This method is called Godunov splitting. If the two subproblems are not 
commutable, the splitting error is O(At), see Chapter 17 of reference [20] , 

There is a 3-stage splitting called Strang splitting, which has a leading error 
term 0(At 2 ), 

i AfP AtK iAfP 

jjn eS > jj* e > w * e_- ^ 
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We will see in the section about diffusion that in reality the errors induced by 
the two splittings are very close. That is because the coefficient of the term 
O(At) is much smaller then the coefficient of 0(At 2 ). Thus, usually Godunov 
splitting is sufficient. 

Because the unsplit problem is dominated by the hyperbolicity, the system 
must be driven mainly by the Boundary Conditions (B.C.) through the first 
subproblem. Thus we discuss the B.C. of the hyperbolic part in the next sub- 



section and present the treatment of B.C. for the parabolic part in Section 3.6 
together with Crank-Nicolson scheme. 



3.2 Initial and boundary conditions 
3.2.1 Initial conditions 

Assume we are interested in the blood flow in an arterial segment (0, L) within 
a time interval (0,T). For an evolutionary problem, a proper I.C. is needed. 
In reality, the information contained in I.C. flows out after a certain interval of 
time, and it will not have influence on the system thereafter. Thus, the I.C. can 
be set arbitrarily, say, U(t = Q,x) = (Aq,0), for convenience. 



3.2.2 Inlet and outlet of the hyperbolic part 



Let us look back to the vector Eq. ( 10 ) again. The two components of this 
system are 

dW-\ dWi 

^1+^(10-0, (15a) 
dWo dWo 

sr + -o. < 15b > 

Since the two eigenvalues have opposite signs, there is exactly one incoming 
characteristic at each end of the computational domain. The incoming char- 
acteristic carries information from outside the domain and thus is essential to 
guarantee the problem to be well-posed. That is to say, the system must be 
supplemented by B.C. in the form 

Wi(0,t) = W 2 (L,t) =g 2 (t),t> 0. 

The outgoing characteristic carries information from inside the domain, 
which can be given by the differential equations. Since Wi t 2 are constants 
along the lines D t Xi^(t) — Ai,2 in time-space plane, we can get W.^ +1 (0) and 
W" +1 (L) by interpolation in the data of the n-time step: 

W? +1 (0) = W 2 n (-A2(0)Ai), W[ l+1 {L) = W?(L- Xi(L)At). (16) 



The characteristics are then transformed to physical variables by relation ([12 
for numerical computation. 
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In reality, we rarely have the explicit expression of incoming characteristics. 
Usually, we want to impose B.C. in physical term A, Q or P. At the inlet, if 
A(t) is given, one can use the relation (111 to deduce: 

W? +1 = W? +1 + 8y^A(t)V2. 

If Q(t) is given, similarly we obtain an approximation 



W? +1 = -W? +1 



A(0)- 

If P(t) is given, from the wall relation ^ simplified with no viscous effect 
(y s = 0), we in fact impose: 



'"- J - "-"- J ^8^{P(t) + pAl /2 ). 



W? +l = W. 



At the outlet, some part of the perturbation of outgoing characteristic W\ may 
be reflected, 

W? +1 = W$ - R t (W? +1 - W?) 

where Rt is the coefficient of reflection. If R t = 0, the B.C. is nonreflecting. 
That means that the outgoing characteristic goes out without leaving any effect 
and that the incoming characteristic is a constant in time. If Rt ^ 0, the 
reflection is due to the resistance in the downstream arteries. 



3.2.3 Conjunction points 

There are many cases when conjunction points need to be considered: when 
there are changes of topology, sharp changes in geometry or mechanic properties. 
Topological change corresponds to the large amount of bifurcations and some 
trifurcations in the arterial network. Sharp changes correspond to the sharp 
variation of the properties of the vessel wall, e.g. sharp increase of stiffness K 
due to stenting or Aq due to aneurysm. As the derivatives of the corresponding 
variables in the source term lead to a singularity or very near a singularity, the 
vessel can be treated as two segments conjuncted together. 

Since all of the conjunctions can be treated with the same method, we con- 
sider a branching point as a sample problem: a main vessel with two daughter 
arteries. At the branching point, there are then six boundary conditions, A™ +1 
and Q£ +1 for the outlet of the parent artery and Q^A^ 1 and Q^ +1 

for the inlets of the two daughter arteries. From the physical point of view, the 
scheme has to preserve the conservation of mass flux 

q; +1 -qT-q£ 1 = ( 17a ) 

and conservation of momentum flux 

i (l n+1 i O n+1 

Ip(%t) 2 + Pp +1 -\p(%) 2 -PI +1 = ^ « = 1.2. (17b) 
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Pp +1 and -P^. +1 are given by the constitutive relation M. Moreover, the outgo- 
ing characteristics in each artery can be determined by the interpolation formula 
(16). is given by the interpolation on the n-th step data of the parent 

vessel and it must be equal to Wi(U™ +1 ), which is given by relation (11 1 



(w 1 ); +1 -w 1 (u£ +1 ) = o. 

The same holds for W% on the two daughter arteries. 



{W 2 



\n+l 



W 2 (U t 



r 1 ) 



1.2. 



(17c) 



(17d) 



Combining Eqs. (17a I, (17b), (17c) and (17d), there are 6 Eqs. with 6 un 



knowns. This nonlinear algebraic system can be readily solved by Newton it- 
erative method or some other nonlinear algebraic solvers with U n as the initial 
guess. 



3.3 MacCormack scheme 

In FD framework, MacCormack method [55] is very suitable for nonlinear hyper- 
bolic systems of conservation laws. It is equivalent to the Lax-Wendroff scheme 
for linear systems. It has the following characteristics: conservative, three-point 
spatial stencil and two time levels (predictor and corrector), second-order accu- 
rate in time and space. 



For the conservative system (13), an approximate solution U* is obtained 
from U n in the first step and then corrected in the second to give the solution 
U n+1 at the next time step t + At. The numerical solution is performed in 
a mesh with N+1 points which defines the spatial resolution Ace = jj, see 
Figure [I] The finite difference equations (at the interior grid points) are then : 



1. predictor step 



r; ( ;" -^(F? +1 -FD + AtS?, i = 2,...N 



2. corrector step 



where F* and S* are evaluated as functions of the predicted solution U* . Note 
that the predictor step applies a forward differencing and the corrector step a 
backward differencing. The order of the two kind of differencing can be reversed. 
The grid points i = 1 and i = N + 1 represent the boundary conditions. 



3.4 Taylor-Galerkin scheme 

In this section, we follow the presentations of Formaggia et al. (TOl [11] and 
Sherwin et al. [37] for the Taylor-Galerkin scheme. Other forms are also possible, 
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L 



X, 



Figure 1: Mesh for FD and FE 



see, for example [H]. The Taylor series for U n+1 up to second order is 

dU n At 2 d 2 U n 



u n+1 = u 1 



At- 



&t 



at 2 



•©(AO- 



Prom Eq. (13), one obtains, 



dU n 



S 7 ' 



dF n 
dx 



(18) 



(19) 



Differentiation both sides with respect to t and exchange of the order of spatial 
and temporal differentiation in the second term gives, 



d 2 U n 
dt 2 



(20) 



where Su = §§ and H = §£. Substituting Eq. (fl9[) into Eq. (po| and then 



both of them into Eq. (18), one gets 



i)U 



At 



(21) 



-Ai(£" + — S%S n ) 



For convenience, we adopt the notation 
F LW (U) = F(U)- 



At 

Y 

At 



H(U)S(U), 



Slw{U) = S{U) + —Su(U)S(U). 

The piecewise linear function space associated with the mesh (Figure [T]) is given 
as, 

V? = {[vh?\v h G C°,v h \ [xi , Xi+x] G C\v h (0) = v h {L) = 0,t = 1..JV}. 

The shape function in this space has the property tpi(xj) = Sij, where 5y is 
Kronecker delta. This is both the trial function space and test function space 
in Galerkin framework. U is approximated by Uh G V®. We further define the 
inner product 



(u, v) 



uvdx. 
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Replace U by Uh in Eq. ( |2l[ ), multiply both sides by basis test functions, and 
integrate over the domain [0, L], finally we get 

(ur\^) = (uj:^ l ) + At(F L MK)^)~(s u (uj:)^^^) ^ 

U^ +1 and are expanded as Uh — Y?j=2 Ujifrj- Instead of evaluated directly 
as nonlinear functions of U£, the terms F(U%), F LW (U^), S LW (U^), Su(U%) 
and H(U 7 h l ) are projected onto the trial function space and expanded by a group 
finite element method. That is, for example, F(UJ^) = X)i'=2 ^JVj> with F™ = 
F([/™). Finally, the matrix form of the FE scheme writes 

At 2 At 2 

MU n+1 = MU n + AtK T F£ w + AtMS n LW - —MF n - —ICF n , (23) 



where 

M i:j = (ipi,ipj), JCij = (A, -q 1 ) 

and 



dx ' dx 



Please note our abuse of notation that in this equation U n+ , F n etc. stand for 
discretized vectors whose values are associated with the grid points. Also note 
that M and JC are functions of S u and H, therefore they must be updated in 
every time step. 

3.5 Second order finite volume scheme 

For finite volume method, the domain is decomposed into finite volumes or cells 
with vertex Xi as the center of cell x i+i/2\> see Figure [2] For every cell, 

the conservation law must hold, 

asi+i/a QJJ rx l+1/2 Qp rx i+1 / 2 

-j—dx+ I -p—dx— I Sdx. 

Application of Gauss's law on the second term gives 

QJJ rx l+1/2 

dx + F\ Xi+x/2 -F\ Xi _ 1/2 = Sdx. (24) 

In the cells, average values are considered, 

Ui = — / U(x)dx, S t = — / S(x)dx. 

Ax Jxi-if* Ax Jxt-1/2 
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Thus Eq. 24 turns into an ordinary differential equation 

(25) 



dt Ax 

We have a local Remiann problem at each interface of neighboring cells, since 
the values of U at the two sides of the interface, Ul to the left and U R to the 
right, are not equal in general. By solving the Remiann problem, a numerical 
flux F* can be obtained. Depending the approximate approaches on solving 
the Remiann problem, different numerical fluxes are possible. Among them, 
Rusanov (or called local Lax-Friedrichs) flux is widely used. It writes 

F(U L ) + F(U R ) U R -U L 
F {U L , U R ) = c , 

with 

c= sup (sup |Aj(?7)|) 

U=U L ,U R j£l,2 

where Xi(U) and \2(U) are the eigenvalues of J c . HLL flux is another option 
and it has less numerical diffusivity. Since Rusanov flux is more simple and 
robust, it is adopted in this paper. HUl and U R equal the average values at the 
cell, the numerical flux will be of first order accuracy. Linear reconstructions 
of U within the cells are necessary for a second order numerical flux. To this 
end, let us consider the slope of a scalar s within the z-th cell Dsi, which can be 
approximated by (si — Sj_i)/Ax, (s^+i — Si)/Ax or (s,+i — Sj_i)/2Ax. Then 
the values of s at the interfaces associated with this cell can be recovered as 

Ax , Ax 

s i-i/2+ = St — -Dsi and s l+1/2 _ = s, + ~Y Ds i- 

The discretization of derivative in space can achieve a second order accuracy 
by this method. But the solution will have nonphysical oscillations. Some ex- 
amples of oscillations induced by these methods can be found in Chapter 6 of 
reference |20) . Slope or flux limiter and non-oscillatory solution are integral 
characteristics of FV schemes. MUSCL (monotonic upwind scheme for conser- 
vation law) is one popular slope limited linear reconstruction technique. To 
present MUSCL, we first define a slope limiter, 

!min(x,y) if x, y > 0, 

max(x,y) if x,y < 0, 

else 

Then the slope Dsi is modified as 

7-, . 1 S &i \ 

Dsi = mmmod : , : . 

y Ax Ax ' 

The values of A and Q at the interface can be obtained as 

Ax Ax 

A l -i/ 2 + = Ai- —DAi, A i+1/2 - =A t + —DAi 
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and 

Air Ax 
Q1-1/2+ = Qi ^DQi, Qj+1/2- = Qi + —DQi. 

It is easy to verify that the variables are conserved by this reconstruction 
Ri-l/2+ + Ri+1/2- _ „ Qi-l/2+ + Qi+1/2- _ n 

The adopted numerical integration in time is also of second order accuracy. Let 



us rewrite Eq. ( 25 1 as 

f -v ■» 

where 

Aa; 

Note that the fluxes have been replaced by numerical ones. A 2-step second 
order Adams-Bashforth (A-B) scheme is applied for the temporal integration, 

c/r +i = u? + At(^(un l^ur 1 ))- 

This scheme can be initiated by a forward Euler method. A second order Runge- 
Kutta (R-K) method is also possible. But the R-K method requires once more 
resolution of $(£/) at every step. This may be offset by a larger time step allowed 
by the R-K method. But note also that the boundary conditions are determined 
dynamically. The A-B method allows less resolutions of the nonlinear algebraic 
equations at conjunctions points. Thus we choose the A-B method for the 
temporal integration. 

x ;-i/2 X i+U2 L 



Figure 2: Mesh for FV 



3.6 Treatment of the parabolic subproblem 

For the previous 3 schemes, only the hyperbolic subproblem resulted from split- 
ting is solved. For the parabolic subproblem, Crank-Nicolson method is very 
suitable. The temporal and spatial discretization has the form, 

ui L+1 - u[ = c v ( u?£ - zu? +1 + u^ 1 u* +1 -zu* + uu ^ 

At 2 [ Ax 2 Ax 2 h 

where U* is the solution of the first hyperbolic subproblem. The matrix of 
the resulting algebraic system is tridiagonal, which is quite cheap to invert. 
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This scheme is second order accurate both on time and space. Moreover, it is 
unconditionally stable. It is natural to set a homogeneous Neumann B.C. for 
the parabolic subproblem, d x U p (0,t) = d x U p (L,t) — 0. The subscript p stands 
for parabolic. 



3.7 Local Discontinuous Galerkin scheme 

In FV framework, the recovery of Ul and Ur of higher accuracy requires bigger 
stencil. In higher dimensions, this kind of reconstruction leads to difficulties. 
On the other hand, it is quite straightforward to increase the order of approxi- 
mation polynomials in one finite element. Unlike the global FE, the neighboring 
elements do not share the same values at the interface. Numerical flux are ob- 
tained from these values, where the dynamics of the system can be considered. 
We present a nodal DG scheme, following Hesthaven and Warburton's book [T5] . 
The domain is decomposed into K non-overlapping elements, see Figure [3} At 
each element, the local approximation to the solution is a polynomial of order 
N = N p — 1. The global approximation to U is direct summation of these local 
solutions: 

k=K 

rk 



(26) 



k=l 



Similarly, the flux F and the source term S can also be approximated by di- 
rect summation of piecewise N-th degree polynomials. The local form of the 
conservation law on the fc-th clement has the form, 



du£ 
dt 



dFt 
dx 



si 



(27) 



Multiplication by a test function 
over one element gives 



at both sides of Eq. ( 27 ) , and integration 



dt ' 



dFt 

dx 



Applying integration by part on the second term, we have: 



(28) 



dul 
dt ■ 



Ft 



dx 



Xk + l 



(29) 



Again, since the value of Ut at the two sides of the interface, U^(xk) and 
Ufr +1 {xk), are not equal, a numerical flux F* is introduced here. Through the 
numerical flux, information is communicated between elements. In practice, the 
second term is integrated by part again for convenience of computation. Thus 
we have 



■dvl 

dt : 



dFt 

dx 



+ (-Ft;cb k + F*<j> K ) 







A 


Xk 







(30) 
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If we introduce N p nodes within the element the local solution can be 



expanded as 



u£(x,t) 



N p 

E 

i=l 



(31) 



where l k (x) is the Lagrange interpolant associated with the z-th point. For the 
Galerkin scheme, Eq. (30) must hold for every test function l k {x). Thus we have 
N p equations for N p unknowns. In matrix form, the system can be written as, 



M 



k dU k 
dt 



K. k F k + (-F k r + F*r) 



M k S k , 



(32) 



where 



(It I 



3 ) D k 



dl k 
t ul i jk\ 



and l k is the vector of equations (l\ , l k , -./^ ) T - The system of equations can be 
turned into a semi-discrete form, 



dt 



(M k Y l {-F k l k + F*l k 



k ik 



Xk+l 



s k . 



(33) 



where 



V 



(id) 



-1 y k 

Hid) 



dr 



is the local differentiation operator [15]. The computation of M k and T> k is 
crucial. We define an affine mapping from a reference element (—1, 1) to D k , 



x(r) — x k + 



1 



(x k +i - Xk). 



The local operators can be readily computed as, 
M k = J k 



kljdr, 



dlj 
Jk—r- 
dr 



where = {x k+ i — x^)/2 and ij is the Lagrange interpolant at the reference 
element. Note that the terms f -^liljdr and ^-| can be precomputed and 
stored. Legendre-Gauss-Lobatto points has to be chosen as the interpolation 
points to minimize the computation error. For more details, we refer to Chapter 
3 of reference [TS]. For the temporal integration, a second order A-B scheme is 



applied for reasons as discussed in Section 3.5 



The scheme previously presented can treat a hyperbolic problem. But in 
this setting Crank-Nicolson method is hard to apply, because the values at the 
interfaces are duplicated. We consider the problem formulation of Eq. ([5]), where 
the flux contains convective part F c and diffusive part F v . For the convective 
part, Godunov flux is applicable. For the diffusive flux, a straight idea is to use 
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the central flux, (F v (Ul) + F v (Ur))/2. But as pointed out by Shu el al. 
this choice is inconsistent. 

To solve this problem, we rewrite the original equations as 



dU , d{F c ~C v q) 



dt 



Ox 



= s 



dQ n 

ox 



In semi-discrete form, the equations for one element are 



dt 



= -T> k F k - (M k y l (-F k l k + F*l k ) 



k\-li 



Xk+l 



<l 



V k Q k - {M k )- l {~Q k l k + Q 



The flux in these equations have to be modified accordingly: F k = F k — C v q k , 
F* = F* — (C v q)*. F* is defined by Godunov flux. The numerical flux (C v q)* 
and Q* are defined by the central flux. The introduction of an auxiliary variable 
q stabilizes the scheme. Note that the auxiliary equation does no involve time 
evolution, and thus q k can be eliminated in every time step. The addition 
of storage or computational cost is very limited. This method is called local 
discontinuous Galerkin scheme. 







D, 



Figure 3: Mesh for DG 



4 Results and discussion 

In this section, the numerical solutions in various cases are verified against ana- 
lytic, semi-analytic solutions or clinical observations. At first, the computations 
are done on a single uniform vessel. In case of small perturbations, a linearized 
system is obtained. If this system is homogeneous, it allows pure wave solution. 
If the source terms due to the skin friction and the viscosity of the wall are added 
respectively, asymptotic solutions are obtained. In case of larger perturbations, 
the full nonlinear system allows shocks. The shock-capturing property of each 
scheme is tested in this case. After the test on a single vessel, a simple bifurca- 
tion is computed and the reflection and transmission coefficients arc compared 
with analytic ones predicted by linearized system. At the final of this section, 
a network with 55 arteries is computed and the numerical solutions are checked 
against clinical observations reported in literature. 
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4.1 Propagation in an uniform tube 

In this subsection, we compare the numerical results with analytic ones for 
a pulse wave on a single uniform vessel (d x ((3\fA^) — d x (3 — 0). To avoid 
reflections, non reflective B.C. is set at the outflow to mimic a semi-infinite 
tube. Adding a small perturbation ((eA, eQ)) to the equilibrium solution (U = 
(j4o,0)), substituting it into the governing equations and dropping the terms 
with quadratics of e, we obtain the equations for the perturbations in a linear 
form: 

d_A 8Q_ dQ d_A _C l ~ d*Q 

dt + dx _0 ' dt +Ca dx~ A Q Q + Cv d x i (34) 



with c = J£^/Aq, the Moens-Korteweg celerity. To investigate the propa- 
gation phenomena at first, we drop the non-homogeneous part (Cf = and 



C v = 0). Then Eqs. (34) become d'Alembert equations, which admit the 
pure wave solution. If we assume that the initial condition is at equilibrium 
(A = A , Q = 0), and the inflow is prescribed as Q(0,t) = Qi n {t) with 

2ir T 
Q in (t) = Q c sin(—t)H(-t + y ), t> 0, 

where H (t) is the Heaviside function, T c the period of the sinusoidal wave and 
Q c the amplitude. The solution is CqA = Q = Qi n {x — cot), which means that 
the waveform propagates to the right with a speed of Cq. 

We propose a numerical test with parameters of the tube inspired by |37j : 
L = 250cm, A = 3.2168cm 2 , /3 = 1.8734 x 10 6 Pa/m and p = 1.050 x 10 3 kg/m 3 , 
and accordingly, cq = 400cm/s. To impose a small perturbation, we choose 
Q c = lml/s and T c = 0.4s. In this case the change ratio of the radius is 
AR/Rq = Qo/(2A c ) = 0.04%, thus the perturbation is assured small enough. 
We take the linearized analytic solution at time t = OAs as reference to compute 
the errors of the numerical solutions. The norm of the error vector is defined by 



\E\ 



1 I Qnumerical Qanalytic 



To specify the time step, we define the formula, At = Ctj^r- For the DG 
scheme, the time step formula is modified accordingly as At = ^fj^r, with V 
the degree of the polynomial. The value of CV must satisfy the CFL ^Courant- 
Friedrichs-Lewy) condition. 

To test the spatial convergence, we fix C t — 0.1, and vary the number of 



mesh nodes N. The log-log plot of \\E\\ against Aa; can be seen in Figure 4(a) 



We have two main observations. First, all of the schemes converge with an order 



between 1 and 2 and the DG schemes converges faster (see Figure 4(a) ). Second, 
as shown by Figure |4(b)| the difference between the analytic solutions and all of 
the numerical solutions are hardly discernible with a moderate number of mesh 
points (N TG = N FV = N FD = 800, N DG - Vl = N DG --p 2 = 100). 
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To test the temporal convergence, we fix the mesh (Ntg = ^fv = ^ra = 
800, Ndg-Vx = Ndg-v 2 = 100) and increase (or decrease) the time steps, 



the error varies slightly for all of the schemes except FV (see Figure 4(c)). If 
Ct > 0.6, Taylor-Galerkin, FV and MacCormack become unstable. For the DG 
schemes, C t can not be greater than 0.1 (see Figure [4(d) \ . For the convergence 
of the temporal integration, the FV scheme has to choose a smaller time step 
than the value prescribed by the CFL condition. 

To compare the actual speed and accuracy of the four schemes, we set iV 
and Ct such that the errors achieve the same order of magnitude (see Table [I] 
and Figure [4(f) \ . Except the Taylor-Galerkin scheme, all of the schemes have 
the similar accuracy with very close running time (see Figure [4(e) and |4(f)[ ). At 



this point, the Taylor-Galerkin scheme shows the worst accuracy and needs to 
run the longest time. We note that large global matrix arise in Taylor-Galerkin 
scheme while the operators in other schemes are local and have small size. That 
explains the relative poor performance of Taylor-Galerkin although a larger 
time step is allowed by this scheme, we will see that in case of a network of real 
size, the largest number of N is about 100 and Taylor-Galerkin shows a good 



balanced property between accuracy and speed (Section 4.6) 



scheme 


N 


c t 


Taylor-Gakerkin 


800 


0.5 


FV 


800 


0.3 


MacComack 


1600 


0.5 


DG-Vx 


200 


0.1 


DG-V 2 


100 


0.1 



Table 1: Number of elements and coefficient of time step 



4.2 Attenuation due to the viscosity of blood 



We now consider the same linearized system Eq. (34) with the small term due 
to skin friction (C/ ^ and C v — 0). The main dynamics of the system 
will be grossly the same traveling wave but attenuated by viscosity of blood. 
This behaviour can be predicted by asymptotic analysis. We have a small non- 
dimensional parameter e/ = T c Cf/Ao, which is the ratio of the characteristic 
time of pulse T c to the characteristic time of attenuation A /Cf. In order to see 
how the waveform slowly evolves when it propagates to, say right, we make a 
change of variables to r = e/f and £ = x — cgt (slow time, moving frame). The 
two differential operators dt and d x expand as 

d dr d 
dt = dtdr 



<9£ d 


d 


d 


dtd£ 




d£ 


d 


d£ d 


d 


dx 


dx <9£ 
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(a) 



(b) 




Figure 4: Test on a uniform tube. Top left, errors with different size of elements 
(cells). Top right, all the solutions for the pulse wave at time 0.4s are overlapped 
and the analytic solution is indicated by cross signs. Middle left and right, with 
a fixed mesh (N TG = N FV = N FD = 800, N DG _ Vl = 200 and N DG --p 2 = 100), 
errors as functions of C t coefficient. Bottom left and right, running time and 
errors for the configuration shown in Table [T] 
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The solution has the asymptotic expansion 



A = A + e f A x + Q = Q + (fQi + ■■■ 

Substituting these into the governing equations expressed in new variables and 
collecting the terms with the same order of e f , one has 

. dA 9Q , ,dA 8A 1 ^dQi. 
dQ . 2 d Ao, , ,dQ a dQx 2 dA 1 Q . 

We take the first order term in ef in the first equation, substitute it in the first 
order term in e/ in the second equation, finally we obtain 

dQ dA Q 

From the terms of the zeroth order in e/, which involve derivative in £ only, the 
solution must have the form Q — coA (t, £) + 4>{t). Substituting it into the 
previous equation implies terms ^ and 4>(t). These are secular terms and thus 

can be set null. So we have coA = Q and = — ^rQo , or 

Qo = Oo(0,0^ r/(2Tc) = Qo^x-co^e'^/^. 

For more on asymptotic analysis of blood flow in large blood vessels, we refer 
to reference [45] . 

In Figure [5] we plot the snapshots of the waveform at time 0.2s, 0.4s, 0.6s 
and 0.8s. In the computation, the initial and boundary conditions are the same 
as in the previous subsection. The mesh and the time steps in Table [l] are 
adopted. The damping rate of the amplitude of the waveform agrees very well 
with the analytic prediction, exp(— 2 Aoc )' wmcn ^ s indicated by the dashed line. 
Also note that the errors of different schemes are not the same. The FV scheme 
causes the peak of the wave to slightly flatten, while all of the other schemes 
are dispersive: we have small oscillations at the foot of the signal. 

4.3 Diffusion due to the viscosity of the arterial wall 



We now consider the same linearized system Eq. (|34j) but with the Kelvin- Voigt 
effect and no viscous fluid effect (Ct = and C v 7^ 0). The small parameter 
is now e v = C v /(cqT c ). If we apply the same technique as described in the 
previous subsection, we can readily obtain the diffusive behaviour of the pulse 
wave in the moving frame: 

dQ _ c\T c d 2 Q 

dr "2 W 1 ' 
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Taylor-Gale rkin 

FV 

MacCormack 

DG 




Figure 5: Attenuation due to the skin friction. The snapshots are at time 0.2s, 
0.4s, 0.6s and 0.8s. The dashed line is exp(-^|;) with 2A c Q /Cf ~ 2000cm 
The flux is normalized with respect to Q c . 



The solution of this equation can be given by the convolution 



Qo(r,0 



?o(0,OG(r^-C)dC 



where G is the fundamental solution of the heat equation (35) 

1 



V^tc 2 T c 



-«7(2rc^T c ) 



and Qo(t, £) is the initial state. In the test vessel, the parameters are kept 
the same as in the case of attenuation. The coefficient C v is 0.6275m 2 /s and 
e v ~ 0.1. This corresponds to 0=5OOOPa • s, which is in the range of observed 
values on animals [2] • To facilitate the calculation of the analytic solution, non- 
reflecting B.C.s are imposed at the two ends of the vessel and I.C. is a half 
sinusoidal waveform for flux Q (dashed line in Figure [6]) and a constant cross 
section Aq . It is clear that half of the initial wave propagates to right and 
at the same time the waveform is spread out due to the diffusive effect. The 
analytic solution at time 0.4s (indicated by cross signs) agrees well with the 
corresponding numerical solution. 

Another point worthy noticing is the operator splitting errors. In the DG 
scheme, no operator splitting error is induced. All of the other numerical 
schemes adopt operator splitting method. They produce very accurate solutions 
as well as DG. Thus it verifies the a priori judgement that Godunov splitting 
is sufficient. 



4.4 Shock-like phenomena due to the nonlinearity 

We now consider the full nonlinear system, but without any source terms (C/ = 
and C v — 0). The small parameter is now ti — Qc/(cqAq). If we apply the 
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Figure 6: Diffusion due to the viscosity of the wall. The dashed line is the 
initial condition. One half of the original waveform propagates to right. The 
snap shots are at time 0.2s, 0.4s, 0.6s and 0.8s. The analytic prediction from 
the convolution at time 0.4s is indicated by cross signs. The difference between 
the different numerical solutions is not discernible. The flux is normalized with 
respect to Q c . 



same technique as described in the previous subsection, we can readily obtain 
the equation for the non linear behaviour of the pulse wave in the moving frame 
(inviscid Burgers equation): 



dQo _ dQo 
dr ~ 2A ^° 3i 



(36) 



One important consequence of nonlinear hyperbolic system is that shocks 
may arise even if the initial condition is very smooth. In normal physiological 
conditions, shocks are not observed in arterial systems. But in venous system, 
shock-like phenomena may occur on muscular veins during walking and running. 
The intramuscular pressure (equivalent to P ext in our model) can rise to 20 — 40 
kPa in a few ms [3]. In such situation, experiment and numerical simulations [HI 
I25j have shown this critical behaviour. For another example, the traumatic 
rupture of the aorta is responsible for a significant percentage of traffic death 
and the rupture may be well accounted for by the shock-like transition resulted 
from the blunt impact to the thorax [TB]. Thus we test all of the schemes 
in such case. To generate a shock, only two parameters are modified: L = 
800cm and Q c = 200ml/s. The change ratio of the radius is about 7.78%. The 
number of elements for Taylor-Galerkin, FV and MacComack schemes is 800. 
The DG scheme uses 200 elements and the order of polynomial is 2. Figure [7] 
shows that a shock starts to appear near the point 300cm. Strong oscillations 
are generated at the front foot of the waveform by Taylor-Galerkin scheme. 
On the other hand, strong oscillations are induced at the back of the large 
gradients. For the DG scheme, there are some smaller oscillations both in front 
and back. That is because the characteristic structures are taken into account 
in the numerical flux. Limiters may be introduced to eliminate the oscillations. 
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This remedy will be necessary for DG to be applicable on problems with shocks. 
For the FV scheme, the shock is well captured without nonphysical oscillations. 
That verifies the total- variation-diminishing (TVD) property of the MUSCL FV 
scheme. 

On Figure |7(b)| we plot a case with some viscosity of the wall. The added 
moderate physical diffusive term smoothes the oscillations and all of the schemes 
give almost the same result. 




Figure 7: Shock due to the nonlinearity of the system. The left figure (a) shows 
that sharp gradient arises in a nonlinear hyperbolic system from smooth initial 
condition. Numerical schemes may cause spurious oscillations. FV scheme with 
a flux limiter captures the shock without nonphysical oscillations. The right 
figure (b) shows that all of the schemes give almost the same result for a system 
with a moderate physical diffusive term. 



4.5 Reflection and transmission at a branching point 

Up to now, we focused on the various behaviours of wave within a single vessel: 
propagation, attenuation, diffusion, etc. Now, we look at the boundaries of each 
artery. Indeed, pressure waves are reflected and transmitted at the conjunction 
points of a network. For a linearized system, given the impedance Z — the 
reflection and transmission coefficients at a branching point can be calculated 
by the formula: 

n = Z p 1 ~ ( Z ch + Z d 2 1 ) j- = 2Z p 1 

Zp 1 + (z£ + Z^) ' Zp- 1 + (Z^ 1 + Z^) ' 

where Z p and Zd are the characteristic impedance of the parent and daughter 
vessels [T31 [3T] . In Figure [8j for sake of illustration, the configuration of the 
branching and the time profile of pressure at two locations are shown. The 
amplitude is normalized with respect to Q c = 1 x 10~ 6 m 3 /s = lml/s. For the 
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parent vessel: j3 — 2.3633 x 10 6 Pa/m, Ao = 4cm 2 and for each of the daughter 
vessels: f3 — 6.3021 x 10 6 Pa/m, Ao — 1.5cm 2 . According to the formula, 
1Z = 0.2603 and T = 1.2603. The pressure profiles at the points A and B 
agree very well with the analytic predictions. All of the numerical schemes are 
compatible with this treatment of conjunction. 
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Figure 8: Reflection and transmission of pressure wave at a branching point. 
The time profiles of the pressure at points A and B are plotted. The analytical 
reflection and transmission coefficients are 0.2603 and 1.2603 (indicated by the 
dashed line). 



4.6 Application on a full systematic arterial network 

As already mentioned in the introduction, a relatively realistic description of 
arterial system has been done in ID simulations, with different numerical solvers 
by different teams. For example, in [551137], Galerkin approach is used. In these 
papers, wall viscosity is not included. Note that [35] gives a survey of literature 
on the details the model, and adopted a viscoelastic model of the wall. But, 
in all of those papers, usually only one numerical scheme is adopted and cross 
comparisons among them are not available. In this subsection, we compute a 
network of 55 arteries with the viscoelastic model we presented above and make 
a cross comparison among the numerical schemes. To this end, the topology 
and properties value of the arterial network are adapted from [37] . But the 
viscosity coefficient of the Kelvin- Voigt model on human body is not given in 
this paper. In reference [2], the viscosity of aortic wall of dogs was modeled 
by a Kelvin- Voigt model and it shows that the value of <j> is in the range of 
3.8 ± 1.3 x 10 3 Pa • s to 7.8 ± 1.1 x 10 3 Pa • s. Hence, we assume 4> = 5 x 10 3 Pa • s 
to calculate the coefficient C v . The final parameters of the network we use are 
shown in Table [3j We note that there may be differences between arteries in 
human and dog and the arteries in different locations may cause a considerable 
variation. Nevertheless the inclusion of viscosity term makes it possible to test 
the numerical schemes in a more realistic condition. 

The peak value of the input flux Q c is 500 ml/s. This value is very close to 
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the peak flow rate at the root of aortic artery [35]. We choose min-f | 5 (L*/cg) 

as a reference length, with L % the length of the i-th artery and Cq the wave 

speed of the linearized system. For a coarsest possible mesh, the number of 

elements (cells) of each artery is N£ ase = L mill i=55( L i/ e i) J > where |_-J is floor 

function. We computed the relative change of solutions when the number of 

the elements (cells) are doubled. Figure [9] shows the relative change of the 

solutions when the number of the elements (cells) is changed from 2Ni, ase to 

4-/Vbase- The relative change of a quantity (for example flux Q) with two meshes 

Ni and N 2 is defined as \ \Qn 1 -Qn 2 \ \\/{Qmax -Qmin), where ||- ||i stands for 1- 

norm, Q ma x and Q m in are the maximum and minimum values within one heart 

beat. Figure [9] shows that the change of flux and pressure are less than 1.5% 

for all of of the schemes except DG. Thus we plotted in Figure [10] the results 

computed with mesh 2Nj )ase . The DG scheme is not tested in this manner 

because it is already converged: result in Figure [T0| shows that the error for DG 

is already very small with the coarsest possible mesh. Time step is prescribed 

by At = Ct min ]z ? 5 ( jk— ) . The coefficient C t and the corresponding real time 

. 0. 1 — 1 
step in the computation is shown in Tabic Y2\ 




Figure 9: Relative change of the solutions when the mesh is doubled from 2Nt> ase 
to ANb ase . The left figure shows that the relative change of all of the flux is less 
than 1.3 %. The right figure shows that the relative change of all of the pressure 
is less then 0.6% . 



In Figure 10 we plot the history profile of flux and pressure at the middle of 
four representative arteries. All of the numerical solutions agree very well. The 
main features of the pressure and flux profiles reported in literature |37[ 135] are 
observed. The peak value of pressure waveform increases as we travel down the 
system. We can also see the dicrotic notch at artery 1. At artery 37, a reverse 
flow is observed (see Figurc [T0(f)[ ), which agrees with clinical measurement [55] . 
The result in this paper is smoother than the result on the corresponding arteries 
in [37]. The result with viscosities is closer to the clinical observations [35] . We 
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scheme 



C t At (10 6 s) running time (min) 



Taylor-Gakerkin 



0.4 222 22.0 

0.25 139 31.9 

0.1 55.5 91.2 

0.006 6.66 576 



FV 
MacComack 
DG 



Table 2: Time steps and running time for one heart beat using one processor 
on a standard Linux work station. 

realize that it is very important to consider the wall viscosities to give more 
realistic result. This agrees with the conclusion drawn by the comparison of 
numerical results with data from in vitro experiments [TJ [35] . 

5 Conclusions 

In this paper, we incorporated a Kelvin- Voigt viscoelastic constitutive relation 
of arterial wall with a ID blood flow model. This led to a hyperbolic-parabolic 
system which was then solved by four numerical schemes: MacCormack, Taylor- 
Galerkin, second order finite volume and discontinuous Galerkin. The imple- 
mentations were verified with analytic, semi-analytic or clinical observations 
in many cases. At first, a single uniform tube was considered. Under the as- 
sumption of small nonlinearities, we obtained asymptotic solutions of linearized 
systems with different source terms. The propagation, attenuation and diffusion 
of the waveform were illustrated by both the numerical and analytic solutions. 
Moreover, in case of a larger nonlinearity, the shock capturing property of each 
scheme was tested. After the test on a single vessel, a simple bifurcation was 
computed to check the numerical coupling of different arteries. Finally, we com- 
puted a relatively realistic network with 55 arteries. The check of the numerical 
solutions in all cases was very favorable for all of the schemes. The schemes 
can be compared in four aspects: the accuracy, the ability to capture shock 
phenomena, the computation speed and the complexity of the implementation. 

1. MacCormack and Taylor-Galerkin schemes generate small oscillations. FV 
scheme has slight arbitrary steepening effect. Both diffusion and disper- 
sion errors are very small for DG. Nevertheless all of the schemes converge 
with a moderate fine mesh and precisely capture the various phenomena 
of this hyperbolicity-dominated hyperbolic-parabolic system. 

2. MacCormack and Taylor-Galerkin perform very poorly when there is a 
steep gradient. Both of them present strong oscillations at one side of the 
jumping location. DG scheme has smaller oscillations at both sides of the 
jump. Numerical flux limiters are possible to filter out the oscillations. 
That will further complicate the schemes and the theory and technique 
is still under research [T71 [35] . On the other hand, there are very mature 
techniques to impose a slope limiter in the FV scheme. Shock capturing 
property is unique for FV among the four schemes presented in this paper. 
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Flux 



Pressure 




(h) (i) 

Figure 10: The history profiles of flux and pressure at four locations. Ten heart 
beats are computed to secure steady state is achieved and the tenth heart beat 
is plotted. The differences between the four numerical schemes are very small. 
See Table [3] for time steps and running time of each scheme. 



3. For a network of human size, the speed of computation can be ordered from 
fast to slow as: Taylor-Galerkin, second order FV, MacCormack and DG. 
The temporal integration in the Taylor-Galerkin scheme is more efficient 
than Adams-Bashforth 2-step method. Thus it allows a larger time step 
with a comparable accuracy. But if the number of elements for one artery 
is too large (larger than 500) , the speed of Taylor-Galerkin becomes slower 
because the size of the global matrix increases quadratically and thus 
the storing and inverting of matrix become expansive. The DG scheme 
prevents the application of Crank-Nicolson method on the diffusive term. 
An explicit method called local DG scheme was adopted in this paper. 
Even with a moderate diffusion coefficient (within the range observed in 
physiological condition), a very small time step is necessary for stability. 
That makes the computation of 1 heart beat cost about 9 hours while all 
of the other schemes cost only 20 to 90 minutes (using one processor on a 
standard Linux work station). 

4. From easiest to hardest, the implementation of the schemes can be ordered: 
MacCormack, second order FV, Taylor-Galerkin and local DG. 

As a final conclusion from the point of view of practical application, we 
recommend MacCormack in case of small nonlincarities as it is very simple 
and robust; second order FV will be a very good option if there maybe shock- 
like phenomena in the systems; Taylor-Galerkin has quite balanced properties 
between speed and accuracy if no shock-like phenomena may present in the 
system; DG is suitable for systems with very small physical diffusive term since 
both the numerical diffusion and dispersion are very small in this scheme. 

Our future work may include these directions: including a physical dispersive 
term to account for the axial tension of vessels; parallelizing the scheme to 
speedup the computation and coupling the arterial system with the venous 
system. 
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Table 3: Arterial network 







I 


A 





Cv 




ID 


Name 


(cm) 


(cm z ) 


(10 6 Pa/cm) 


(10 4 cm 2 /s) 


Rt 


1 


Ascending aorta 


4.0 


6.789 


0.023 


0.352 


— 


2 


Aortic arch I 


2.0 


5.011 


0.024 


0.317 


- 


3 


Brachiocephalic 


3.4 


1.535 


0.049 


0.363 


— 


4 


R. subclavian I 


3.4 


0.919 


0.069 


0.393 


— 


5 


R. carotid 


17.7 


0.703 


0.085 


0.423 


— 


6 


R. vertebral 


14.8 


0.181 


0.470 


0.595 


0.906 


7 


R. subclavian II 


42.2 


0.833 


0.076 


0.413 


— 


8 


R. radius 


23.5 


0.423 


0.192 


0.372 


0.82 


9 


R. ulnar I 


6.7 


0.648 


0.134 


0.322 


— 


10 


R. interosseous 


7.9 


0.118 


0.895 


0.458 


0.956 


11 


R. ulnar II 


17.1 


0.589 


0.148 


0.337 


0.893 


12 


R. int. carotid 


17.6 


0.458 


0.186 


0.374 


0.784 


13 


R. ext. carotid 


17.7 


0.458 


0.173 


0.349 


0.79 


14 


Aortic arch II 


3.9 


4.486 


0.024 


0.306 


— 


15 


L. carotid 


20.8 


0.536 


0.111 


0.484 


— 


16 


L. int. carotid 


17.6 


0.350 


0.243 


0.428 


0.784 


17 


L. ext. carotid 


17.7 


0.350 


0.227 


0.399 


0.791 


18 


Thoracic aorta I 


5.2 


3.941 


0.026 


0.312 


— 


19 


L. subclavian I 


3.4 


0.706 


0.088 


0.442 


— 


20 


L. vertebral 


14.8 


0.129 


0.657 


0.704 


0.906 


21 


L. subclavian II 


42.2 


0.650 


0.097 


0.467 


— 


22 


L. radius 


23.5 


0.330 


0.247 


0.421 


0.821 


23 


L. ulnar I 


6.7 


0.505 


0.172 


0.364 


— 


24 


L. interosseous 


7.9 


0.093 


1.139 


0.517 


0.956 


25 


L. ulnar II 


17.1 


0.461 


0.189 


0.381 


0.893 


26 


intercoastals 


8.0 


0.316 


0.147 


0.491 


0.627 


27 


Thoracic aorta II 


10.4 


3.604 


0.026 


0.296 


— 


28 


Abdominal aorta I 


5.3 


2.659 


0.032 


0.311 


— 


29 


Celiac I 


2.0 


1.086 


0.056 


0.346 


— 


30 


Celiac II 


1.0 


0.126 


0.481 


1.016 


— 


31 


Hepatic 


6.6 


0.659 


0.070 


0.340 


0.925 


32 


Gastric 


7.1 


0.442 


0.096 


0.381 


0.921 


33 


Splenic 


6.3 


0.468 


0.109 


0.444 


0.93 


34 


Sup. mensenteric 


5.9 


0.782 


0.083 


0.439 


0.934 


35 


Abdominal aorta II 


1.0 


2.233 


0.034 


0.301 


— 


36 


L. renal 


.3.2 


0.385 


0.130 


0.481 


0.861 


37 


Abdominal aorta III 


1.0 


1.981 


0.038 


0.320 


— 


38 


R. renal 


3.2 


0.385 


0.130 


0.481 


0.861 


39 


Abdominal aorta IV 


10.6 


1.389 


0.051 


0.358 


— 


40 


Inf. mesenteric 


5.0 


0.118 


0.344 


0.704 


0.918 


11 


Abdominal aorta V 


1.0 


1.251 


0.049 


0.327 




42 


R. com. iliac 


5.9 


0.694 


0.082 


0.405 




43 


L. com. iliac 


5.8 


0.694 


0.082 


0.405 




44 


L. ext. iliac 


14.4 


0.730 


0.137 


0.349 




45 


L. int. iliac 


5.0 


0.285 


0.531 


0.422 


0.925 


46 


L. femoral 


44.3 


0.409 


0.231 


0.440 




47 


L. deep femoral 


12.6 


0.398 


0.223 


0.419 


0.885 


48 


L. post, tibial 


32.1 


0.444 


0.383 


0.380 


0.724 


49 


L. ant. tibial 


34.3 


0.123 


1.197 


0.625 


0.716 


50 


L. ext. iliac 


14.5 


0.730 


0.137 


0.349 




51 


R. int. iliac 


5.0 


0.285 


0.531 


0.422 


0.925 


52 


R. femoral 


44.4 


0.409 


0.231 


0.440 




53 


R. deep femoral 


12.7 


0.398 


0.223 


0.419 


0.888 


54 


R. post, tibial 


32.2 


0.442 


0.385 


0.381 


0.724 


55 


R. ant. tibial 


34.4 


0.122 


1.210 


0.628 


0.716 



Data adapted from |37| and [2], 
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